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ABSTRACT 

Precision measurement of the scalar perturbation spectral index, n s , from the Wilkinson Microwave 
Anisotropy Probe temperature angular power spectrum requires the subtraction of unresolved point 
source power. Here we reconsider this issue, attempting to resolve inconsistencies found in the lit- 
erature. First, we note a peculiarity in the WMAP temperature likelihood's response to the source 
correction: Cosmological parameters do not respond to increased source errors. An alternative and 
more direct method for treating this error term acts more sensibly, and also shifts n s by ~ 0.3er 
closer to unity. Second, we re-examine the source fit used to correct the power spectrum. This fit 
depends strongly on the galactic cut and the weighting of the map, indicating that either the source 
population or masking procedure is not isotropic. Jackknife tests appear inconsistent, causing us 
to assign large uncertainties to account for possible systematics. Third, we note that the WMAP 
team's spectrum was computed with two different weighting schemes: uniform weights transition to 
inverse noise variance weights at I — 500. The fit depends on such weighting schemes, so different 
corrections apply to each multipole range. For the Kp2 mask used in cosmological analysis, we prefer 
source corrections A = 0.012 ± 0.005 /iK 2 for uniform weighting and A = 0.015 ± 0.005 /iK 2 for iV bs 
weighting. Correcting WMAP's spectrum correspondingly, we compute cosmological parameters with 
our alternative likelihood, finding n s = 0.970 ± 0.017 and as = 0.778 ± 0.045. This n s is only 1.8a 
from unity, compared to the ~ 2.6a WMAP 3-year result. Finally, an anomalous feature in the source 
spectrum at I < 200 remains, most strongly associated with W-band. 

Subject headings: cosmology: observations, cosmic microwave background, cosmological parameters, 
methods: data analysis 



1. INTRODUCTION 

Measuring n s , the spectral index of initial scalar fluctu- 
ations, which is scale invariant (n s = 1) in the Harrison- 
Zeldovich model and slightly shallower in inflation mod- 
els, is difficult, primarily because experimental systemat- 
ics require control over a broad range of spatial scales. In 
inflation, the deviation from unity closely relates to the 
inflationary potential and the number of e-folds of ex- 
pansion, so a statistically robust measurement of n s =/= 1 
places compelling constraints on the physics of the infla- 
tionary epoch. 

Because all-sky measurements of the cosmic microwave 
background (CMB) access the largest observable scales 
in the universe, the angular power spectrum of the CMB, 
with a long lever arm, is crucial to such studies. Indeed, 
the latest data release from the Wilkinson Microwave 
Anisotropy Probe ( WMAP) claims ~ 2.6a deviation from 
the Harrison- Zeldovich spectrum (jSpergel et al.l [20071 ). 
Unfortunately, the CMB is not a totally clean measure- 
ment. For example, the well-known degeneracy with 
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the optical depth since recombination (t) makes preci- 
sion measurement of n s impossible using CMB temper- 
ature anisotropics alone, and polarization is required to 
break it. Complicated noise properties and hints of un- 
known systematics in the WMAP measurement of large- 
scale polarization indicate that the systematic uncer- 
tainty in both t and n s sh ould still be considered sig- 
nificant (jEriksen et al.ll2007n . 

Another important, but under-appreciated, complica- 
tion for the measurement of n s is additional power in the 
angular spectrum from unresolved, and unmasked, point 
sources. At high this shot noise can significantly bias 
the power spectrum, and consequently n s . The WMAP 
team devised a sensible prescription for dealing with this 
contaminant: 1) Use the spectral energy distribution 
measured from detected sources (and distinct from the 
CMB) to infer it for undetected ones; 2) measure the 
contamination using multifrcquency data; 3) correct the 
spectrum; and 4) marginalize over the measurement er- 
ror w hen computing the likelihood (|Hinshaw et al.ll2003l . 
l2007h. 

iHuffenberger et ail (|2004l ) found a level of source con- 
tamination c onsistent with the lev el in the first WMAP 
data release (|Hinshaw et al.ll2003[). However, based o n 
the three year temperatu re data (jHinshaw et al.l f2007h . 
IHuffenberger et all (|2006f ) measured a point source spec- 
trum with two irregularities. First, at / > 200 the spec- 
trum is white, but with an amplitude b e low t he value 
in the original preprint of iHinshaw et alJ ()2007l ). In the 
present work, we di scovered a small er r or in the power 
spectra used for the IHu ffenberge r et all (|2006l ) estimate, 
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which should have reported A = 0.013 ± 0.001 /jK 2 in- 
stead of A = 0.011 ± 0.001 //K 2 , still below the original 
WMAP va lue of A = 0.017 ± 0. 002 /zK 2 . Prompted by 
our result, iHinshaw et al.l |2007) rc-cxamincd the issue, 
revising their value down somewhat and increasing the 
error b ars, to A = 0.014 ± 0.003 ^K 2 . The lSpergel et all 
(|2007f ) bispectrum analysis indicates a non-Gaussianity 
consistent with these values, but lacks the statistical 
power of the multifrequency power spectrum comparison. 
T he second peculiarity is tha t the power at 100 < I < 200 
in Huffcn berger et all (|2006f ) was inconsistent, at strong 
statistical significance, with the rest of the white spec- 
trum. 

This paper again considers the power spectrum source 
correction procedure in detail. We begin in Section [2] 
with a study of the impact of the source correction on 
the scalar spectral index through the likelihood. Follow- 
ing this, we probe the source amplitude in Section [3l ex- 
amining the dependence of the fit on the sky weighting, 
mask, year of observation, and frequency dependence, 
and present our best estimates of the cosmological pa- 
rameters. These same tests probe the robustness of the 
I < 200 feature. Finally, wc conclude in Section [U 

2. SOURCE CORRECTION IMPACT ON SPECTRAL INDEX 

The final WMAP temperature power spectrum is a 
noise-weighted combination of cross spectra computed 
from V- (61 GHz) and W-band (94 GHz) maps. Prior 
to the computation of the angular spectrum, a fore- 
ground model is removed from the maps. The individual 
cross-spectra are corrected for the sky mask, instrument 
beams, a nd point source cont amination before combina- 
tion (see IHinshaw et aT1l2007h . The combined spectrum 
is folded into the likelihood calculation, 8 which interfaces 
to a Markov Chain M onte Carlo code such as CosmoMC 
(jLewis fc Bridldl2002[ ). yielding parameter estimates. 

Assuming the angular spectrum for sources is white, 
we wish to explore the dependence of n s on the size of 
the source corre ction. In Figu r e [H w e show the cor- 
rection made bv IHinshaw et al.l (|2007l ). The correction 
is not white, because the combined spectrum gets more 
contribution from W-band at higher /, due to the noise 
weighting. To change the amplitude of the correction, 
we simply scale. Later, we discuss the estimation of this 
amplitude. 

We explore two methods for incorporating the source 
errors into the likelihood. In the first, we use the stan- 
dard WMAP likelihood method. The total log-likelihood 
is split in two, —2 logi = C = £q+£i, where Co is due to 
cosmic variance, noise, and mask-induced mode coupling, 
and Ci includes t he uncertainty due to the p oint source 
and beam errors (IHinshaw et al.ll200l l2007f) . The com- 
putation of C\ assumes that total likelihood is Gaussian, 
and uses the Woodbury formula to compute an update to 
the likelihood. Under the Gaussian assumption this split 
is exact. However, the £o term is not treated as Gaus- 
sian, but comp uted using the Ga ussian plus log-normal 
approximation (jVerde et al.ll2003l ). 

Treating the source and beam uncertainty separately is 
very fast (since the uncertainty can be well approximated 
by a small number of modes), but not really necessary. 

8 We are using version 2.2.2 of the WMAP likelihood, available 
at http://lambda.gsfc.nasa.gov/. 
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Fig . 1. — The size of the source correction in IHinshaw et aLl 

It avoids the inversion of a ~ 1000 x 1000 matrix per like- 
lihood evaluation, but the low-^ part of the code already 
inverts a ~ 2000 x 2000 matrix. Performing this inver- 
sion and including the beam and source term with the 
other sources of error in Cq is therefore little additional 
burden. 

We take this step in our second method for incorpo- 
rating the source errors. We have modified the WMAP 
likelihood code to integrate the beam/point source co- 
variance matrix into the cosmic variance/noise/mask co- 
variance matrix. (This requires inverting the original 
Fisher matrix, adding the beam/point source term, and 
inverting back.) Then we feed the new covariance matrix 
to the Gaussian plus log-normal approximation. Under 
this procedure, the change in the likelihood due to the 
inclusion of the beam/point sources term is -2.64, com- 
pared to C\ = —1.22 computed with the Woodbury ex- 
pansion, for the test theory spectrum included with the 
WMAP likelihood code (where C ~ 3541). 

In the top panel of Figure[2l we show the dependence of 
n s on the source correcti on. We hold the errors fixed at 
the IHinshaw et all (|2007f) value of a A = 0.003 [lY 2 : , and 
marginalize over all the other parameters. As the source 
correction increases, the power spectrum at high I is low- 
ered, and n s decreases. Completely ignoring the source 
correction shifts n s by ~ 0.01 higher, or about 0.6<t. Our 
alternative likelihood procedure shifts n s higher by 0.005, 
or 0.3cr. 

Next we study the influence of error in the source cor- 
rection. At very small and at very large error, we expect 
the parameter measurement to be independent of the 
error. For very small source error, the errors in other 
quantities dominate. For large errors, all modes which 
could be contaminated by point sources are effectively 
projected out, and the parameter measurement is again 
independent. Near S/N ~ 1, the measurement should 
undergo a transition, where the measurement error in- 
creases and (possibly) the mean value changes. In the 
bottom panel of Figure [21 we show the dependence of 
n s on the error on the source correction amplitude. The 
WMAP split likelihood method shows an unexpected re- 
sult. The size of the point source error apparently has 
no effect on the measurement of n s . Even with the point 
source error rivaling the size of the acoustic peaks, there 
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Fig. 2. — (Top) The change in n B for given source corrections, 
A. We compare the WMAP team's likelihood code to our mod- 
ified version, which computes the high-i likelihood in a slightly 
different way. The e rror in the s ource correction is fi xed at the 
IHinshaw etaTJ H2007fl value used in lSpergel et al.l 020071). The thin 
lines bound the 68% probability interval. The lSpergel et al.l II2007T ) 
values are slightly offset horizontally for visibility, and either ignore 
(higher n s ) or include (lower) a correction for Sunyaev-Zeldovich 
effect contamination. The other n s values ignore SZ. (Bottom) 
Using the WMAP team's source correction, this plot shows the 
dependence of n s on the point source error, ua- 

is no effect. The values of the likelihood change, but not 
the distribution of points in the Markov Chain. This 
seems to be a clear indication that there is something 
wrong. The dependence of n s on A in the top panel im- 
plies that at least the error bars should increase as a a 
increases. 

On the other hand, our modified likelihood shows the 
expected behavior. As the source correction error is 
made very large, the errors increase by about 38 percent, 
and the mean value moves above 1. For this likelihood, 
the modes not subject to contamination by sources actu- 
ally prefer n s > 1, and a solid measurement of the source 
contamination is vital to the measurement of n s < 1. We 
also note that the errors on the source measurement do 
not make much difference, as long as a a < 0.003. In 
the next Section, we examine this measurement in more 
detail. 

3. UNRESOLVED POINT SOURCE SPECTRA 
3.1. Method 



The point source spectrum can be estimated via a lin- 
ear combination of the individual cross-spectra at several 
frequencies, a combination which projects ou t the CMB 
component. iHuffenberger et aTl (|200 4. 2006) examined 
the unresolved source component in WMAP data, using 
a ge neralized version of t he method for the same task 
from IHinshaw et all (|2003f ). 

In addition to V- and W-band used for cosmological 
measurement, the source analysis uses Q-band (41 GHz) 
because the contaminating sources are much brighter at 
lower frequencies. There are 276 cross-spectra in these 
three bands, accounting for all combinations of differenc- 
ing assemblies per channel, and treating the three years 
of observation separately. These are combined as 

(1) 



to give the point source amplitude in a band de- 
noted by L, where the weight is W% for the binned 
cross spectrum estimate D 1 ^ , made from maps i and 

j. The map cross spectra D % [ include the CMB power 
spectrum and the contribution from sources. By virtue 
of being cross-spectra (i ^ j), they are noisy but do 
not have a noise bias. The weights are based upon the 
frequency independence of the CMB signal (in thermo- 
dynamic temperature units), the spectral energy distri- 
bution of the sources (measured for bright sources as 



-2.0, S oc yP +2 (IHinshaw et al.l l200l iTrushkinl 



l2003t IHinshaw et al.l l2007f) ). and the estimated noise 
covariance in the cross s pectru m measurements. (See 
IHuffenberger et all (|2004l I2006D for details.) 
The weights obey the constraint 



0, (2) 

L' ij 

which means that no CMB will leak into the point source 
estimate if the maps are properly calibrated and the 
instrumental beams arc perfectly deconvolved from the 
spectra. If the source spectral energy distribution is cor- 
rect, the weights also provide an unbiased estimate of the 
point source power spectrum even if the noise covariance 
is wrong (although an incorrect noise covariance leads to 
sub-optimal estimates and incorrect error bars). 
The covariance of the source power spectrum estimate 

is 

(A L A L ,) = (WE D W) LLS (3) 

where Sd is the estimate of the cross spectrum covari- 
ance matrix. If the covariance matrix is diagonal in cross- 
spectra and multipole bin, then the weights are diagonal 
in multipole. Under this assumption, we plot some ex- 
ample weights in Figure [3] 

The expected shot noise angular spectrum of sources 
in the WMAP data is flat in Cj, so we plot our measured 
spectra as I versus C;. Throughout we normalize the 
spectrum at Q band in antenna temperature, which gets 
a larger signal than V or W. 

In the following subsections, we describe several tests 
of the point source measurement to explore the robust- 
ness of these features and their origin. Tests included 
changing the weights on the map for computing the spec- 
trum, modifying the mask, fitting for the galactic hemi- 
spheres separately, fitting year by year, and changing the 
assumed source spectral energy distribution. 
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Fig. 3. — The weights for various cross-spectra in the source 
measurement, using the Kp2 mask, and two different map weight- 
ing. The individual weights for the 276 cross-spectra were summed 
together based on frequency coverage, yielding the 6 combinations 
shown. 

3.2. Spectrum at I > 200 

Whatev er the cause of the excess p ower at I < 200 re- 
ported bv iHuffenberger et al.1 (|2006f ) may be, it is very 
likely not relevant for correcting the spectrum at high 
I. In this Section, wc therefore first concentrate on the 
correction required for the cosmological analysis, consid- 
ering only I > 200, and then return to the anomalous 
low-/ feature in Section EOl 

In a maximum likelihood estimate of the power spec- 
trum, the map is weighted by the pixel-pixel inverse sig- 
nal plus noise covariance. At smal l scales this estimate is 
computationally impractical, and iHinshaw et alj (|2007f ) 
instead approximate it using two weighting schemes in 
the computation of the power spectra. In the signal 
dominated regime at I < 500, they use spectra where 
every pixel is weighted evenly; in the noise dominated 
regime at I > 500, the maps are inverse noise weighted 
(i.e. weighted by the number of observations, iV b s ). 

We compute the cross-spectrum using these two 
schemes and estimate the source contribution, plotted 
in Figure 01 The flat portion of the spectrum is notably 
higher in the N Q b s weighted case. Specifically, for a white 
noise fit including I > 200, A = 0.015 ± 0.001 /xK 2 ver- 
sus A = 0.012 ± 0.001 /jK 2 for the flat weighting. (The 
single amplitude fits for all combinations we tried are 
compiled in Table Q]) This leads immediately to a key 
point: Since the combined spectrum is built from two 
weighting schemes at different Z's, the source correction 
must be different as well. Thus we should have a smaller 
source correction for I < 500 and a large r correction for 
I > 50 0. The weighting scheme used by IHinshaw et al.l 
(|2007t ) is not specified. 

This source level difference is a strong indication that 
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Fig. 4. — Fits for the unmasked source component, comparing 
spectra computed from maps using uniform weight outside of the 
mask with using N \, a (inverse variance) weighting, plotted as the 
Q-band amplitude. The lowest /-bin has no detection and is not 
shown. 



TABLE 1 
Fits for source amplitude 



Mask Weight Subset A a a X 2 /dof 

(ixK 2 , ant. Q) 



Diagonal Covariance 



30kp0 


fiat 


all 


0.006 


0.001 


0.25 


kp2 


fiat 


VW 


0.006 


0.004 


1.99 


kpO 


fiat 


all 


0.008 


0.001 


0.54 


kp2 


flat 


noQf 


0.010 


0.001 


2.02 


kp2 


fiat 


noQ2 


0.010 


0.001 


0.71 


kp2 


flat 


yrl2 


0.011 


0.002 


1.02 


kp2 


fiat 


yrl3 


0.011 


0.002 


1.10 


kp2 


fiat 


all (N. hemi.) 


0.011 


0.001 


0.92 


kp2 


fiat 


all 


0.012 


0.001 


1.39 


kp2 


flat 


yr22 


0.012 


0.003 


1.16 


kp2 


flat 


yr23 


0.012 


0.002 


0.50 


kp2 


fiat 


QW 


0.013 


0.001 


0.70 


kp2 


fiat 


all (S. hemi.) 


0.013 


0.001 


0.68 


kp2 


fiat 


noWlW2 


0.013 


0.001 


2.48 


kp2 


fiat 


noW3W4 


0.013 


0.001 


1.23 


kp2 


fiat 


QV 


0.014 


0.001 


1.65 


kp2 


nobs 


all 


0.015 


0.001 


0.69 


kp2 


fiat 


yrll 


0.017 


0.003 


0.87 


kp2 


fiat 


yr33 


0.017 


0.003 


0.20 


Diagonal + beams errors 


kp2 


fiat 


noQf 


0.009 


0.001 


1.87 


kp2 


fiat 


VW 


0.010 


0.004 


2.07 


kp2 


fiat 


noQ2 


0.010 


0.001 


0.73 


kp2 


fiat 


yrl2 


0.011 


0.002 


1.05 


kp2 


flat 


yrl3 


0.011 


0.002 


1.10 


kp2 


fiat 


yr22 


0.011 


0.003 


1.11 


kp2 


fiat 


yr23 


0.011 


0.002 


0.34 


kp2 


fiat 


all 


0.012 


0.001 


1.53 


kp2 


flat 


QW 


0.013 


0.001 


0.72 


kp2 


nobs 


all 


0.015 


0.001 


0.79 


kp2 


flat 


QV 


0.016 


0.001 


2.30 


kp2 


flat 


yr33 


0.017 


0.003 


0.25 


kp2 


flat 


yrll 


0.019 


0.003 


0.85 



Note. — Fits for source power spectrum amplitude, 
including data at I > 200. 
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Fig. 5. — The difference map, Q— fLC, shows two point-like 
objects near the LMC, left, compared to the Q-band noise rms 
for the same region, right, based on the number of observations 
reported by the WMAP team. 

the source population contaminating the WMAP spec- 
trum is not isotropic. One possible scenario is the fol- 
lowing: Because source positions are stochastic, and 
source power strongly favors brighter sources (C; oc 
/ dS S 2 dN/dS), brighter sources which, by chance, fell 
in the best observed regions could boost the spectrum 
in the iV b s weighting over uniform weighting. How- 
ever, the size of this effect is much too small to be 
a viable explanation. Wc computed power spectra of 
Monte Carlo realizations of a noiseless map containing 
isotropically distributed faint sourc es, based on differen- 
tial source counts dN/ dS oc S~ 2 3 ([White fc Majumdau 
12004 IClearv et al.ll2005h . normalized to WMAP source 
counts at 1 Jy, setting an upper flux limit to reproduce 
a reasonable amount of power. Over 1000 simulations, 
for a bin 100 < I < 150, the rms fluctuation of the differ- 
ence in power between the uniform and iV b s spectra is 
1.8 percent, and should be smaller at higher I. Thus the 
odds are very slight that the change in the power is due 
to chance alignments of sources with the well observed 
part of the sky; sources at the appropriate flux level are 
too numerous. Alternatively, the observed anisotropy in 
the source population could be cither be representative of 
the real sources, or indicate a problem with the masking 
procedure, with bright sources slipping through. 

The iVobs weighting emphasizes the ecliptic poles, so 
these are the best places to look for suspect sources. In 
particular, visual inspection of a Q— ILC difference map 
(Figure [5]) shows two bright sources near the Large Mag- 
ellanic Cloud, in the highly observed portion of the sky 
near the south ecliptic pole. These sources are as bright 
as some of the sources found by the WMAP source- 
detection algorithm, and indeed t he brighter one is in- 
cluded in t he subsequent c atalog s of lLopez-Caniego et alJ 
(12001 and lNie fc Zhanel (pOOl . based on WMAP data. 
The presence of the nearby LMC may have interfered 
with the source finding procedure. 

Changing the sky cut gives further indications of an 
anisotropic source population. In addition to the Kp2 
cut used in the cosmological analysis, we recomputed the 
point source fit using cross-spectra generated with two 
additional masks and uniform weighting: the KpO mask 
and a conservative mask consisting of the union of KpO 
with a \b\ < 30° cut galactic cut (see Figure [6]). As seen 
in Table [TJ the source power drops significantly as we 
expand the masks, indicating that unmasked sources are 
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Fig. 6. — (Top) Three masks used in this analysis. The smallest, 
and most aggressive mask is the WMAP Kp2 mask, intermediate 
is KpO, and most conservative is the union of KpO and a |6| < 
30° cut. (Bottom) The source power spectrum, computed with 
three different masks. As the mask becomes larger, the power goes 
down, suggesting a concentration of unmasked source power near 
the galactic plane. 

brighter or more common near the plane. This raises two 
possibilities. Either some of the sources are galactic in 
origin, or sources are less efficiently found and masked 
near the plane. 

From a visual inspection of the Q— ILC difference map, 
wc find ~ 6 more bright point-like objects, many near the 
galactic plane, and all in the southern hemisphere (see 
Table [5]). Many of these objects are also already noted 
in the the lLopez-Caniego et alJ (|2007t ) and lNie fc Zhanel 
(|2007t ) catalogs. Computing the power spectrum in hemi- 
spheres (with the Kp2 mask), the source amplitudes 
straddle the value for the whole sky. At I < 200, the 
Northern hemisphere has slightly less source power and 
the the Southern slightly more (See Figure [7]). 

Th e WMAP source mask, bui lt from a variety of cata- 
logs ([Hinshaw et al.ll2003l 120071 ). may not cover the sky 
evenly, because of differences in sensitivity in different re- 
gions. It may be more difficult to identify sources amid 
large galactic foregrounds. This would be particularly 
true for flat spectrum sources in lower frequency radio 
surveys. The WMAP source mask has noticeable gaps 
near the galactic plane, but outside the Kp2 and KpO 
masks. 

Two jackknife tests, breaking the data into subgroups, 
also show some peculiarities. First, we divided up the 
cross spectra based on the years of observation (Fig- 
ure [5]). Each individual year accounts for 28 cross- 
spectra, while a pair of years accounts for 64. At I > 200, 
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Fig. 9. — The source power spectrum estimate for combinations 
of bands. 



TABLE 2 
Sources found in Q— ILC 



RA 


Dec. 


Gal. Ion. 


Gal. lat. 


ID 


00 43 14 


-73 16 30 


303.7 


-44.0 


PMNJ0047-7308V 


01 04 11 


-72 08 02 


301.5 


-45.0 


PMNJ0059-7210*,t 


04 51 56 


-69 31 10 


276.2 


-33.7 




05 20 56 


-66 10 24 


281.0 


-35.6 


PMNJ0506-6109*,+ 


06 49 57 


-16 56 12 


227.9 


-8.0 


PMNJ0650-1637* 


20 50 59 


+28 51 35 


72.5 


-9.7 


_t 


20 52 02 


+31 55 18 


75.1 


-8.0 




20 57 54 


+31 29 08 


75.6 


-9.2 


t 



Note. — Point-like objects, found visually in a Q— ILC differ- 
ence map. Marked sources were fou nd respectively in the 



catalog of lLopez-Caniego et all l|2007f ) and lNie fc Zhangj (120071 ) 
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Fig. 8. — Year by year comparisons of the point source power 
spectrum. 

the fits for individual years are always as greater than 
for pairs of different years. The largest estimate is the 
(year 3 x year 3) fit, at A = 0.017 ± 0.003 /zK 2 , while 
the smallest is (year 1 x year 2), at 0.011 ± 0.002 /zK 2 . 
This seems unlikely to be due to chance, and could have a 



number of causes. For example, a slight cross-correlation 
between Ql and Q2 or between Q and V, introducing a 
small noise bias in the cross spectra, could have this ef- 
fect. Removing cither Ql or Q2 from the source estimate, 
the amplitude drops to A = 0.010 ± 0.001 ^K 2 . 

However, this dependence on the Q band could also be 
due to the shape of the source spectrum. For the second 
jackknife test, we broke the cross spectra up by band. We 
considered combinations with Q and V bands only, Q and 
W only, and V and W only. These spectra are shown in 
Figure [51 At I > 200, the V+W combinations are very 
noisy, and have the lowest source fit, A = 0.006 ± 0.004 
/zK 2 , though the \ 2 ' 1S poor. The Q+V combination has 
the largest source amplitude, A = 0.016 ± 0.001 /jK 2 . 
This discrepancy could mean that the source spectrum 
is incorrect: If the true spectrum is steeper than the 
[3 = —2 we have used, or the spectrum steepens between 
V and W, one could observe this effect. In this case, the 
Q+V combination gives the most accurate measurement 
of the Q band source contamination, but this is not used 
in the cosmological analysis. To get the correction in V 
and W bands, we are better off to scale from the lower 
amplitude, although even this is not completely satis- 
factory. Using the wrong spectrum for sources means 
that the shape of the source correction will be some- 
what incorrect. We gradually steepened the source spec- 
trum, and re-fit. The Q+V, Q+W, and V+W fits can 
be brought within lcr by setting (3 = —2.5. These can 
be made equal (at A = 0.012 fiK 2 ) by setting (3 ~ -2.9, 
which compared to f3 = —2.0 would represent a 50 per- 
cent smaller amplitude in V and an 80 percent smaller 
amplitude in W. If true, the source correction should be 
made correspondingly smaller. (We continue to quote 
values of A in Q-band assuming (3 — —2.0 scaling to V 
and W.) 

In the appendix, we discuss the impact of beam errors. 
These tend to have little impact on the spectrum at I > 
200, since the resulting CMB leakage is large where the 
CMB is large, at lower I. 

Finally, we must decide by how much to correct the 
combined WMAP spectrum. There is a strong detection 
of unresolved power in the WMAP spectra, but as we 
have seen through the above tests, our knowledge of the 
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TABLE 3 

COSMOLOGICAL PARAMETERS 



Parameter WMAP Source/likelihood corrected 



WMAP data only 



n h h 2 


0.0223 ± 0.0007 


0.0220 ± 0.0008 




0.237 ± 0.034 


0.242 ±0.037 


h 


0.735 ± 0.032 


0.731 ± 0.034 


T 


0.088 ± 0.030 


0.090 ± 0.030 


ns 


0.951 ±0.016 


0.968 ±0.017 


CT8 


0.742 ± 0.051 


0.780 ± 0.052 




WMAP + ACBAR H 


- BOOMERANG 


n h h 2 


0.0232 ±0.0007 


0.0224 ± 0.0007 




0.233 ± 0.034 


0.234 ± 0.032 


h 


0.739 ± 0.033 


0.742 ± 0.032 


T 


0.088 ± 0.032 


0.092 ± 0.030 


n B 


0.951 ±0.016 


0.970 ±0.017 


CT8 


0.739 ± 0.051 


0.778 ± 0.045 



Note. — Comparison of marginalized paramet er re- 
sults obtained from Table 5 of lSpergel et al.l (|2007l ) (sec- 
ond column) and after our modifications to the source 
correction and likelihood (third column). 



character of these sources is poor. Our best fit ampli- 
tudes including all the data are A = 0.012 /iK 2 for flat 
weighting and A = 0.015 /iK 2 for -/V b s weighting. But 
the spectral effects (from considering V and W without 
Q) indicate it might be much lower, while individual year 
fits are much higher. We then are forced to the unsat- 
isfying course of artificially inflating the error bars be- 
yond their nominal statistical values, in order to account 
for these possibilities. Therefore we set the error on the 
source estimates at a a = 0.005 /iK 2 . 

To obtain our final estimate for cosmological param- 
eters, we correct the spectrum for point sources us- 
ing the two amplitudes quoted above (uniform weights 
for I < 500 and iV b s weights for I > 500), and our 
modified likelihood code. The marginalized parameters 
from the resulting Markov chains are given in Table 
[3l both considerin g WMAP data alone and including 
data from ACBArIKuo et al.ll2004fl and BOOMERanG 
flJones et a l. 2006; Mo ntrov et al.|[2006t iPiacentini et al.l 
|2006[) . Because of the new likelihood and the broadened 
error bars, this new value of n s = 0.970 ± 0.017 is only 
1.8(7 away from 1. Another consequence is that as in- 
creases to 0.778 ± 0.045. 

3.3. Excess at I < 200 

We now turn to the I < 200 feature, which is incon- 
sistent with a white noise spectrum. This is present in 
both the uniform and A^bs weighted spectra, and in each 
year and pair of years. Because of its shape and promi- 
nence at low I, we initially considered two specific expla- 
nations. The first possibility concerned mis-estimation 
of the overall multiplicative map calibration of each DA 
map: If two maps are calibrated differently, the weights 
in Equation Q] would not cancel the CMB compo- 
nent precisely. Thus, one would observe a leakage from 
the CMB signal into the point source spectrum, with a 
signature resembling the CMB power spectrum. A sim- 
ilar effect would be caused by beam uncertainties. In 
Appendix A, we present the formalism to take these un- 
certainties into account in the method described in Sec- 



tion [331 an d the results from the corresponding analysis 
of the WMAP data arc presented in the bottom Section 
of Table [1] The conclusion from these computations is 
that neither calibration nor beam uncertainties can ex- 
plain this effect. 

A third hypothesis is residual galactic foregrounds, 
which should show through the mask and frequency de- 
pendency. As the galactic cut widens from Kp2 to the 
wide cut (Figure 0, the bin from / = 100-150 drops 
about 0.017 /iK 2 . At the same time the fit for the white 
noise level drops by 0.006 /iK 2 . Subtracting off the white 
noise levels for each, the component in the excess has 
dropped by about 40 percent. The power in the excess 
is still significant, even for this broad cut. 

Next, we turn to the estimates using only two bands. 
In this case, we see that the feature is strongly enhanced 
in the V+W combinations, about the same in the Q±W 
combination, but clearly diminished in the Q+V com- 
bination. The latter appears consistent with the flat 
source spectrum. This may indicate that the excess is 
associated with the W band. On the other hand, be- 
cause the W band spectra tend to carry negative weight 
in these estimates (Figure |3|) , this excess would repre- 
sent a deficit of power in W, which is peculiar. It may 
indicate an over-subtraction of the galactic foreground 
template in W band, which could have consequences for 
the cosmological analysis. At this point, it is difficult to 
be definitive. 

4. CONCLUSIONS 

We have reanalyzed the source correction procedure for 
the three-year data release of WMAP. First, we consid- 
ered the impact of this procedure in the WMAP likeli- 
hood code. Surprisingly, we found that the WMAP likeli- 
hood does not react to changes in the point source correc- 
tion error. We have devised a modified likelihood which 
does respond as expected, although we note that more 
work is needed to completely validate this approach, as it 
couples to the important problem of how to approximate 
a non-Gaussian likelihood with fitting formulas over a 
wide multipolc range. To conclude n s < 1, a precision 
measurement of the source contamination is required. 
We note that the modes not contaminated favor n s con- 
sistent with unity. 

Second, we found several indications that the un- 
masked source population in WMAP data is anisotropic. 
This implies that the combined spectrum should be cor- 
rected differently in two multipole regions, based on the 
weighting of the map. Anisotropy in the unmasked 
sources is unexpected, but can be turned to an advan- 
tage: By very carefully masking near the ecliptic poles 
and galaxy, or employing a wide galactic cut, the point 
source contamination can be cut substantially. This gain 
must be weighed against the reduction of the sky area. 

We note irregularities in jackknife tests of the source 
fit, grouped by time of observation or frequency band. 
This prompted us to adopt large errors on the source fit, 
to account for systematics beyond the estimate of sta- 
tistical error which accompanies our measurement. This 
step is necessary to treat the source correction conserva- 
tively. With the modified likelihood and reduced source 
amplitude from ignoring the excess at I < 200, these en- 
larged error bars are responsible for raising our values of 
n s and a%. 



s 



Finally, the previously noted anomalous I < 200 fea- 
ture is still present, shows signs of being spatially asso- 
ciated with the galaxy, and is most strongly associated 
with the W band. On the other hand, it does not appear 
to be associated with calibration or beam errors. It may 
represent an over-subtraction of the foreground template 
in W, although further investigation is warranted. How- 
ever, the immediate conclusion is that this part of the 
spectrum should not be used to infer the point source 
amplitude at higher Vs. 
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APPENDIX 

LEAKAGE FROM CALIBRATION AND BEAM ERRORS 

In the appendix, we estimate the leakage from CMB into the point source estimate. We consider two sources of 
leakage, the overall calibration of the map and the errors in the beam deconvolution. Each of these can cause the CMB 
spectrum to have a slightly different amplitude or shape in each of the cross spectra. This means that the CMB will 
not cancel itself in the weighted linear combination to give the source estimate. This leakage can be accounted for in 
the cross spectrum covariance. The calibration is a constant function of I, and the beam error is nearly so over a wide 
range of scales, so the shape of the leakage term strongly follows the CMB power spectrum, which drops rapidly with 
I. 

The calibration uncertainty is simpler, so we begin there. In the signal (s) dominated regime, without including 
point sources, we can model the calibration errors with 

d l = (l + g l )s. (Al) 

where the calibration in map di is quantified with dimensionless factors gi with ((/,) = and (<?i<?j) = &g8ij- Then the 
estimated cross-spectra are 

= (\+g l +g J +g l g )C l . (A2) 

Where Ci is the power spectrum of the particular sky realization (although in the actual calculation, we substitute the 
WMAP best- fit spectrum). In the cross spectrum i ^ j, so this is an unbiased estimator, (D¥) = Ci. The covariance 
is 

Cov(£>f , Df, q ) = Cid v [a 2 g (5 ip + 8 iq + S jp + 5 jq ) + a*(5 ip S jq + S iq S jp )] . (A3) 

Essentially, the variance only gets a contribution in terms where the two cross-spectra share a map. 
The point source estimate is 

A i = E H W U' D V =E^'E^^ + * + *9i) (A4) 

V (ij) V (ij) 

where one term has dropped out because the weights sum to zero. Because we have not included any source contribu- 
tion, we have (Ai) = 0, where the ensemble average is over the calibration factors. For plotting, its useful to compute 
the variance due to the calibration factors 

(A?) = J2 C V C V , K'K" [Vg&P + + &3V + + *2(Mi9 + ( A5 ) 

VI" (ij) (pq) 
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Fig. 10. — (Left) Beam errors taken from Ijarosik et al. (2003). The Q-band errors rise smoothly to ~ 0.08 at I = 800. (Right) rms 
leakage of the CMB into the point source measurement, with and without taking the contribution to the beam covariance into account. 

A word about binning the spectrum is appropriate here. In practice we bin the spectrum because the large number 
of cross-spectra (276) and multipoles (~ 800) slows the computation, and the signal-to-noise per multipole in low. We 
may define a binning matrix Gli which averages quantities in non-overlapping bands indexed by L, and compute our 
estimate from the binned cross spectra Y2i GliD^ 3 . The quantity we are now computing is the variance in a bin. Since 
the weights are also computed in bands, this reduces to replacing the Is with Ls and replacing each power spectrum 
term C; with y] ; GuCi- 

In I Jarosik et alj (|2007f ). the calibration uncertainty is quoted as 0.5 percent, so a g = 0.005. As written, we are 
assuming th e all maps are uncorrelated, the case if the calibration fluctuations are dominated by noise as expected 
(|Pagell2007h . In this case, fluctuations in the calibration are fairly small compared to the source estimate. If the 
calibration is dominated by something else, say foregrounds, then whole bands could be correlated, and the calibration 
leakage can be substantial compared to the source fit. 

In a similar way, we can estimate the beam uncertainties. We begin with 

4 m =Bi(l + Ei)s lm , (A6) 

where B\ is an unbiased measurement of the beam, and E\ is the (small) fractional difference between the true beam 
and the measured beam. We assume all beams have independent errors, and define a 5-function-like object 8fp which 
is unity when maps i and j share a beam, and zero otherwise. With this definition, we may use a describe the beam 
error with a Gaussian distribution with (£^) = and 

(EjEj) = >:;/:'e (A7) 
That is, the beam errors are correlated in I but not between beams. The beam-deconvolved cross spectra are 

D\ j = (1 + E})(1 + E>)6 1 . (A8) 
Following the same procedure as with gain calibrations, we find a similar expression for the covariance: 

CovODf , Eft ) = C'i [tf/(6* + O + ^(Sf p + 8f q ) + ^/^(6f p Sf q + Sf q Sf p )] Q, (A9) 

WMAP found that the b eam uncertainty can be well represented as a small number of orthogonal modes (~ 10) 
iHinshaw et all (l200l . [200l : 

^f/ = J2 U rl U rl'- ( Al °) 

r 

where r denotes the mode. Here the binning is a little more complicated than for the calibration uncertainty. Instead 
of simply Ci, the quantities which must be binned are CiU\ (for the first-order terms) and CiU^U^ (for the second 
order term). 

At this time, the beam modes for each WMAP differencing assembly are not public. (Only the modes for the 
comb ined spectrum are in cluded in the likelihood code.) To approximate the beam modes, we take the beam errors 
from Uarosik et all (|2007f ) (shown in Figure [TO]). We treat each of these as a single mode for the associated beam, 
correlating all the multipoles in that beam. This should provide a conservative estimate of the beam covariance, though 
it does not capture all of its properties. With these approximate modes, we can compute the beam variance and the 
rms CMB leakage for any set of weights. This is plotted in the right panel of Figure [TO] For our standard diagonal 
covariance, the rms beam leakage at first looks like a promising explanation for the I < 200 excess. It has a similar 
shape, and strong correlations bin-to-bin. However, when we recompute the weights taking the beam covariance into 
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account, we produce a similar estimate for point sources in a combination which allows very little CMB leakage. Given 
the estimated size of the rms residual CMB, it is surprising that the source estimates are so similar. Perhaps this is 
an indication that the beam covariance is overestimated. In any case, the I < 200 is either not due to beams, or the 
approximation for the beam modes is very poor. 



